<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head><script type="text/javascript" src="https://cdn.jsdelivr.net/gh/opencobra/cobratoolbox@ffa0229fc0c01c9236bb7e961f65712443277719/latest/_static/js/iframeResizer.contentWindow.min.js"></script><meta http-equiv="Content-Type" content="text/html; charset=utf-8"><meta http-equiv="X-UA-Compatible" content="IE=edge,IE=9,chrome=1"><meta name="generator" content="MATLAB 2021a"><title>Convert a reconstruction into a flux balance analysis model </title><style type="text/css">.rtcContent { padding: 30px; } .S0 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 28.8px; min-height: 0px; white-space: normal; color: rgb(213, 80, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 24px; font-weight: normal; text-align: left;  }
.S1 { margin: 20px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: normal; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.S2 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S3 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: normal; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.CodeBlock { background-color: #F7F7F7; margin: 10px 0 10px 0;}
.S4 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 4px 4px 0px 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S5 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S6 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S7 { margin: 10px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S8 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S9 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S10 { color: rgb(64, 64, 64); padding: 10px 0px 6px 17px; background: rgb(255, 255, 255) none repeat scroll 0% 0% / auto padding-box border-box; font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; overflow-x: hidden; line-height: 17.234px;  }
/* Styling that is common to warnings and errors is in diagnosticOutput.css */.embeddedOutputsErrorElement {    min-height: 18px;    max-height: 250px;    overflow: auto;}
.embeddedOutputsErrorElement.inlineElement {}
.embeddedOutputsErrorElement.rightPaneElement {}
/* Styling that is common to warnings and errors is in diagnosticOutput.css */.embeddedOutputsWarningElement{    min-height: 18px;    max-height: 250px;    overflow: auto;}
.embeddedOutputsWarningElement.inlineElement {}
.embeddedOutputsWarningElement.rightPaneElement {}
/* Copyright 2015-2019 The MathWorks, Inc. *//* In this file, styles are not scoped to rtcContainer since they could be in the Dojo Tooltip */.diagnosticMessage-wrapper {    font-family: Menlo, Monaco, Consolas, "Courier New", monospace;    font-size: 12px;}
.diagnosticMessage-wrapper.diagnosticMessage-warningType {    color: rgb(255,100,0);}
.diagnosticMessage-wrapper.diagnosticMessage-warningType a {    color: rgb(255,100,0);    text-decoration: underline;}
.diagnosticMessage-wrapper.diagnosticMessage-errorType {    color: rgb(230,0,0);}
.diagnosticMessage-wrapper.diagnosticMessage-errorType a {    color: rgb(230,0,0);    text-decoration: underline;}
.diagnosticMessage-wrapper .diagnosticMessage-messagePart,.diagnosticMessage-wrapper .diagnosticMessage-causePart {    white-space: normal;}
.diagnosticMessage-wrapper .diagnosticMessage-stackPart {    white-space: normal;}
.embeddedOutputsTextElement,.embeddedOutputsVariableStringElement {    white-space: normal;    word-wrap:  initial;    min-height: 18px;    max-height: 250px;    overflow: auto;}
.textElement,.rtcDataTipElement .textElement {    padding-top: 3px;}
.embeddedOutputsTextElement.inlineElement,.embeddedOutputsVariableStringElement.inlineElement {}
.inlineElement .textElement {}
.embeddedOutputsTextElement.rightPaneElement,.embeddedOutputsVariableStringElement.rightPaneElement {    min-height: 16px;}
.rightPaneElement .textElement {    padding-top: 2px;    padding-left: 9px;}
.S11 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S12 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }</style></head><body><div class = rtcContent><h1  class = 'S0'><span style=' font-weight: bold;'>Convert a reconstruction into a flux balance analysis model </span></h1><h2  class = 'S1'><span style=' font-weight: bold;'>Author: Ronan Fleming, Ines Thiele, University of Luxembourg</span></h2><h2  class = 'S1'><span style=' font-weight: bold;'>Reviewers: </span></h2><h2  class = 'S1'><span>INTRODUCTION</span></h2><div  class = 'S2'><span>Even with quality control during the reconstruction process, it is not appropriate to assume that any reconstruction can be converted directly into a model and used to make predictions. A model must satisfy certain assumptions before it can be used to make reliable predictions. Depending on the type of model model, these assumptions will be different. Each assumption should be chemically or biologically motivated and expressed in an unambiguous manner and preferably both intuitively and mathematically. Flux balance analysis is a mathematical method widely used for studying genome-scale biochemical network. Here one aims to predict steady-state reaction fluxes, where there is a balance between production and consumption of each molecular species that is not exchanged across the specified boundary of a system. In this situation, one might obtain erroneous predictions if the system boundary is incorrectly specified. If a reconstruction contains one or more supposedly mass balanced reactions, but which are actually not mass balanced, such reactions in a model can lead to inadvertent leakage of a metabolite from the model, in violation of mass balance. Similarly, when generating a model for flux balance analysis, it is important to ensure that the network is flux consistent, that is, each reaction can carry a non-zero steady state flux. </span></div><div  class = 'S2'><span>Given a reconstruction with </span><span texencoding="$\hat{m}$" style="vertical-align:-5px"><img src="" width="14" height="19" /></span><span> reactants involved in </span><span texencoding="$\hat{n}$
" style="vertical-align:-5px"><img src="" width="11.5" height="19" /></span><span> reactions, this tutorial demonstrates a method to identify and extract the largest subset of the reconstruction whose internal reactions are both stoichoimetrically and flux consistent and whose external reactions are flux consistent. This model is then mathematically consistent with the basic requirements for generation of predictions using flux balance analysis. The identification of the component of the reconstruction that does not satisfy the aforementioned modelling conditions is also useful for targeting reconstruction effort towards resolving stoichiometric inconsistency or resolving flux inconsistency. The example used in this tutorial illustrates the process of extracting a model consistent with flux balance analsis, from a ReconX reconstruction.</span></div><h2  class = 'S1'><span>PROCEDURE</span></h2><h2  class = 'S3'><span>Select reconstruction to convert into a model and enter parameters</span></h2><div  class = 'S2'><span>Load the ReconX reconstruction, and save the original reconstruction in the workspace, unless it is already loaded into the workspace. </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >clear </span><span style="color: rgb(170, 4, 249);">model</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">if </span><span >~exist(</span><span style="color: rgb(170, 4, 249);">'modelOrig'</span><span >,</span><span style="color: rgb(170, 4, 249);">'var'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(2, 128, 9);">%select your own model, or use Recon2.0model instead</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">if </span><span >1</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        filename=</span><span style="color: rgb(170, 4, 249);">'Recon3D_301.mat'</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        load(filename);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        model=Recon3D;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">else</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        filename=</span><span style="color: rgb(170, 4, 249);">'Recon2.0model.mat'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">if </span><span >exist(</span><span style="color: rgb(170, 4, 249);">'Recon2.0model.mat'</span><span >,</span><span style="color: rgb(170, 4, 249);">'file'</span><span >)==2</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            model = readCbModel(filename);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    model.csense(1:size(model.S,1),1)=</span><span style="color: rgb(170, 4, 249);">'E'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    modelOrig = model;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">else</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    model=modelOrig;</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div></div><div  class = 'S7'><span>Set the level of pri</span><span>nting, zero for silent, higher for more output.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >printLevel=2;</span></span></div></div></div><div  class = 'S7'><span>Choose the directory to place the results</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >basePath=</span><span style="color: rgb(170, 4, 249);">'~/work/sbgCloud/'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%resultsPath=[basePath '/programReconstruction/projects/recon2models/results/reconXs/' model.modelID];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >resultsPath=[basePath </span><span style="color: rgb(170, 4, 249);">'/courses/2019_Leiden_COBRA/practicalsDemo/Day4/' </span><span >model.modelID];</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >resultsFileName=[resultsPath filesep model.modelID];</span></span></div></div></div><div  class = 'S7'><span>Create and enter the folder for the results if it does not already exist</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">if </span><span >~exist(resultsPath,</span><span style="color: rgb(170, 4, 249);">'dir'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    mkdir(resultsPath)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >cd(resultsPath)</span></span></div></div></div><div  class = 'S7'><span>Opt</span><span>ionally create a diary to save the output in case it is very long, this makes it easier to search, especially when debugging the process during the early stages.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">if </span><span >0</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    diary([resultsFileName </span><span style="color: rgb(170, 4, 249);">'_diary.txt'</span><span >])</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div></div><h2  class = 'S3'><span>Ove</span><span>rview some of the key properties of the reconstruction</span></h2><div  class = 'S2'><span>Noting the initial size of the reconstruction is useful for comparisons later with subsets derived according to mathematical specifications.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >[nMet,nRxn]=size(model.S);</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S9'><span style="white-space: normal"><span >fprintf(</span><span style="color: rgb(170, 4, 249);">'%6s\t%6s\n'</span><span >,</span><span style="color: rgb(170, 4, 249);">'#mets'</span><span >,</span><span style="color: rgb(170, 4, 249);">'#rxns'</span><span >)</span></span></div><div  class = 'S10'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="554C34AF" data-testid="output_0" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"> #mets	 #rxns</div></div></div></div><div class="inlineWrapper outputs"><div  class = 'S11'><span style="white-space: normal"><span >fprintf(</span><span style="color: rgb(170, 4, 249);">'%6u\t%6u\t%s\n'</span><span >,nMet,nRxn,</span><span style="color: rgb(170, 4, 249);">' totals.'</span><span >)</span></span></div><div  class = 'S10'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="729932D2" data-testid="output_1" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">  8399	 13290	 totals.</div></div></div></div></div><div  class = 'S7'><span>Make sure the stoichiometric matrix i</span><span>s stored in a sparse format as this accelerates computations with large networks</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >model.S=sparse(model.S);</span></span></div></div></div><h2  class = 'S1'><span>Check in case the reconstruction is a model that is already ready for flux balance analysis</span></h2><div  class = 'S2'><span>There is no need to run </span><span>this live script any further if the reconstruction already satisfies the conditions necessary for flux balance analysis. That is if all internal reactants and reactions are stoichiometrically consistent, and all reactions are flux consistent, then the reconstruction satisfies the criteria to designate it a model ready for flux balance analysis.</span></div><div  class = 'S2'><span>SIntMetBool                     m x 1 Boolean of metabolites heuristically though to be involved in mass balanced reactions.</span></div><div  class = 'S2'><span>SIntRxnBool                     n x 1 Boolean of reactions heuristically though to be mass balanced.</span></div><div  class = 'S2'><span>SConsistentMetBool        m x 1 Boolean vector indicating consistent mets</span></div><div  class = 'S2'><span>SConsistentRxnBool        n x 1 Boolean vector indicating consistent rxns</span></div><div  class = 'S2'><span>fluxConsistentMetBool     m x 1 Boolean vector indicating flux consistent mets</span></div><div  class = 'S2'><span>fluxConsistentRxnBool     n x 1 Boolean vector indicating flux consistent rxns</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">if </span><span >all(isfield(model,{</span><span style="color: rgb(170, 4, 249);">'SIntMetBool'</span><span >,</span><span style="color: rgb(170, 4, 249);">'SIntRxnBool'</span><span >,</span><span style="color: rgb(170, 4, 249);">'SConsistentMetBool'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(170, 4, 249);">'SConsistentRxnBool'</span><span >,</span><span style="color: rgb(170, 4, 249);">'fluxConsistentMetBool'</span><span >,</span><span style="color: rgb(170, 4, 249);">'fluxConsistentRxnBool'</span><span >}))</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">if </span><span >all(model.SIntMetBool &amp; model.SConsistentMetBool)</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            &amp;&amp; nnz(model.SIntRxnBool &amp; model.SConsistentRxnBool)==nnz(model.SIntRxnBool)</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            &amp;&amp; all(model.fluxConsistentMetBool)</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            &amp;&amp; all(model.fluxConsistentRxnBool)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        fullyStoichAndFluxConsistent=1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        fprintf(</span><span style="color: rgb(170, 4, 249);">'%s\n'</span><span >,</span><span style="color: rgb(170, 4, 249);">'Reconstruction is a model that is already ready for flux balance analysis'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">return</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">else</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fullyStoichAndFluxConsistent=0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%s\n'</span><span >,</span><span style="color: rgb(170, 4, 249);">'Reconstruction must be tested to check if it is ready for flux balance analysis'</span><span >)</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div><div  class = 'S10'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="8B3691DA" data-testid="output_2" data-width="420" data-height="18" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Reconstruction must be tested to check if it is ready for flux balance analysis</div></div></div></div></div><h2  class = 'S1'><span>Manually remove certain reactions from the reconstruction </span></h2><div  class = 'S2'><span>Bef</span><span>ore attempting to algorithmically remove stoichiometrically or flux inconsistent supposed internal reactions from a reconstruction to generate a model, there is an option to review the content of the reconstruction and manually identify reactions for removal. That is, there are two options:</span></div><div  class = 'S2'><span>A. Skip manual review of reconstruction content. Move to the next step.</span></div><div  class = 'S2'><span>B. Review the content of the reconstruction and omit any reactions that are assumed to be stoichiometrically or flux inconsistent. With respect to stoichiometric inconsistency, such reactions may be obviously mass imbalanced and not satisfy the heuristic conditions for indentification as an exernal reaction. Alternatively, such reactions may be identified by a previous pass through of this tutorial as being of unknown stoichometric consistent (model.unknownSConsistencyRxnBool(j)==1), after the largest stoichiometrically consistent subset of the network has been is identified. This is an iterative process where multiple rounds of identification of the largest stoichiometrically consistent set and manual curation of the remainder that is of unknown stoichiometric consistency is necessary.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">if </span><span >strcmp(filename,</span><span style="color: rgb(170, 4, 249);">'Recon3.0model'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    modelOrig=model;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">if </span><span >0</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">if</span><span style="color: rgb(14, 0, 255);"> </span><span >1</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span><span style="color: rgb(2, 128, 9);">%Rename some of the biomass reactions to make them more obviously exchange</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span><span style="color: rgb(2, 128, 9);">%reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            model.rxns{strcmp(model.rxns,</span><span style="color: rgb(170, 4, 249);">'biomass_reaction'</span><span >)}= </span><span style="color: rgb(170, 4, 249);">'EX_biomass_reaction'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            model.rxns{strcmp(model.rxns,</span><span style="color: rgb(170, 4, 249);">'biomass_maintenance'</span><span >)}= </span><span style="color: rgb(170, 4, 249);">'EX_biomass_maintenance'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            model.rxns{strcmp(model.rxns,</span><span style="color: rgb(170, 4, 249);">'biomass_maintenance_noTrTr'</span><span >)}= </span><span style="color: rgb(170, 4, 249);">'EX_biomass_maintenance_noTrTr'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span><span style="color: rgb(2, 128, 9);">%ATP hydrolysis is not imbalanced like all the other demand reactions so</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span><span style="color: rgb(2, 128, 9);">%give it a different accronym ATPM = ATP Maintenance</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            bool=strcmp(</span><span style="color: rgb(170, 4, 249);">'DM_atp_c_'</span><span >,model.rxns);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            model.rxns{bool}=</span><span style="color: rgb(170, 4, 249);">'ATPM'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        [model,removeMetBool,removeRxnBool] = manuallyAdaptRecon3(model,printLevel);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">else</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        [model,removeMetBool,removeRxnBool] = manuallyAdaptRecon3Ines(model,printLevel);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    [nMet0,nRxn0]=size(modelOrig.S);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    [nMet,nRxn]=size(model.S);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">if </span><span >nMet0==nMet &amp;&amp; nRxn0==nRxn &amp;&amp; printLevel&gt;0</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        fprintf(</span><span style="color: rgb(170, 4, 249);">'%s\n'</span><span >,</span><span style="color: rgb(170, 4, 249);">'--- Manually removing rows and columns of the stoichiometric matrix----'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        fprintf(</span><span style="color: rgb(170, 4, 249);">'%6s\t%6s\n'</span><span >,</span><span style="color: rgb(170, 4, 249);">'#mets'</span><span >,</span><span style="color: rgb(170, 4, 249);">'#rxns'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        fprintf(</span><span style="color: rgb(170, 4, 249);">'%6u\t%6u\t%s\n'</span><span >,nMet0,nRxn0,</span><span style="color: rgb(170, 4, 249);">' totals.'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        fprintf(</span><span style="color: rgb(170, 4, 249);">'%6u\t%6u\t%s\n'</span><span >,nMet0-nMet,nRxn0-nRxn,</span><span style="color: rgb(170, 4, 249);">' manually removed.'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        fprintf(</span><span style="color: rgb(170, 4, 249);">'%6u\t%6u\t%s\n'</span><span >,nMet,nRxn,</span><span style="color: rgb(170, 4, 249);">' remaining.'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div></div><h2  class = 'S1'><span>Remove </span><span>any trivial rows and columns of the stoichiometric matrix</span></h2><div  class = 'S2'><span>Rem</span><span>ove any zero rows or columns of the stoichiometric matrix</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >modelOrig=model;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >model=removeTrivialStoichiometry(model);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >[nMet0,nRxn0]=size(modelOrig.S);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >[nMet,nRxn]=size(model.S);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">if </span><span >nMet0==nMet &amp;&amp; nRxn0==nRxn &amp;&amp; printLevel&gt;0</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%s\n'</span><span >,</span><span style="color: rgb(170, 4, 249);">'---Checking for Remove any trivial rows and columns of the stoichiometric matrix----'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%6s\t%6s\n'</span><span >,</span><span style="color: rgb(170, 4, 249);">'#mets'</span><span >,</span><span style="color: rgb(170, 4, 249);">'#rxns'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%6u\t%6u\t%s\n'</span><span >,nMet0,nRxn0,</span><span style="color: rgb(170, 4, 249);">' totals.'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%6u\t%6u\t%s\n'</span><span >,nMet0-nMet,nRxn0-nRxn,</span><span style="color: rgb(170, 4, 249);">' duplicates removed.'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%6u\t%6u\t%s\n'</span><span >,nMet,nRxn,</span><span style="color: rgb(170, 4, 249);">' remaining.'</span><span >)</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div><div  class = 'S10'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="00B93836" data-testid="output_3" data-width="420" data-height="18" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">---Checking for Remove any trivial rows and columns of the stoichiometric matrix----</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="85496CEC" data-testid="output_4" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"> #mets	 #rxns</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="A5B173ED" data-testid="output_5" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">  8399	 13290	 totals.</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="644A54A2" data-testid="output_6" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">     0	     0	 duplicates removed.</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="072C7261" data-testid="output_7" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">  8399	 13290	 remaining.</div></div></div></div></div><div  class = 'S2'><span>Che</span><span>ck for duplicate columns by detecting the columns of the  S matrix that are identical upto scalar multiplication.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >modelOrig=model;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >dupDetectMethod=</span><span style="color: rgb(170, 4, 249);">'FR'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >dupDetectMethod=</span><span style="color: rgb(170, 4, 249);">'S'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >removeFlag=0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >[modelOut,removedRxnInd, keptRxnInd] = checkDuplicateRxn(model,dupDetectMethod,removeFlag,printLevel-2);</span></span></div></div></div><div  class = 'S7'><span>Remove any du</span><span>plicate reac</span><span>tions, and</span><span> uniquely involved reactants, from the stoichiometric matrix</span><span>. </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">if </span><span >length(removedRxnInd)&gt;0</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    irrevFlag=0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    metFlag=1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(2, 128, 9);">%set all reactions reversible that are duplicates</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    model.lb(removedRxnInd)=-model.ub(removedRxnInd);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(2, 128, 9);">%remove duplicates</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    model = removeRxns(model,model.rxns(removedRxnInd),irrevFlag,metFlag);</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div></div><div  class = 'S7'><span>Display the statistics on the duplicate reactions, </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >[nMet0,nRxn0]=size(modelOrig.S);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >[nMet,nRxn]=size(model.S);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">if </span><span >nMet0==nMet &amp;&amp; nRxn0==nRxn &amp;&amp; printLevel&gt;0</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%s\n'</span><span >,</span><span style="color: rgb(170, 4, 249);">'---Remove any duplicate reactions----'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    [nMet0,nRxn0]=size(modelOrig.S);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    [nMet,nRxn]=size(model.S);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%6s\t%6s\n'</span><span >,</span><span style="color: rgb(170, 4, 249);">'#mets'</span><span >,</span><span style="color: rgb(170, 4, 249);">'#rxns'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%6u\t%6u\t%s\n'</span><span >,nMet0,nRxn0,</span><span style="color: rgb(170, 4, 249);">' totals.'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%6u\t%6u\t%s\n'</span><span >,nMet0-nMet,nRxn0-nRxn,</span><span style="color: rgb(170, 4, 249);">' duplicates removed.'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%6u\t%6u\t%s\n'</span><span >,nMet,nRxn,</span><span style="color: rgb(170, 4, 249);">' remaining.'</span><span >)</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div><div  class = 'S10'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="D9AD4176" data-testid="output_8" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">---Remove any duplicate reactions----</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="2DAFCBF2" data-testid="output_9" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"> #mets	 #rxns</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="C497C1D8" data-testid="output_10" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">  8399	 13290	 totals.</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="A70DEDD7" data-testid="output_11" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">     0	     0	 duplicates removed.</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="82CCAB47" data-testid="output_12" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">  8399	 13290	 remaining.</div></div></div></div></div><h2  class = 'S3'><span>Rem</span><span>ove any duplicate reactions upto protons</span></h2><div  class = 'S2'><span>Remove reactions reactions that differ only in the number of protons involved as substrates or products. Also remove exclusively involved reactants.</span></div><div  class = 'S2'><span>Save a temporary model for testing, before making any changes.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >modelH=model;</span></span></div></div></div><div  class = 'S7'><span>Find the prot</span><span>on indicies in different compartments. A proton, with index i, is asumed to be represented by an abbreviation within model.mets{i} like h[*], where * denotes the compartment symbol.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >nMetChars=zeros(length(modelH.mets),1);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">for </span><span >m=1:length(modelH.mets)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    nMetChars(m,1)=length(modelH.mets{m});</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >protonMetBool=strncmp(modelH.mets,</span><span style="color: rgb(170, 4, 249);">'h'</span><span >,1) &amp; nMetChars==length(</span><span style="color: rgb(170, 4, 249);">'h[*]'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">if </span><span >printLevel&gt;2</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    disp(modelH.mets(protonMetBool))</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div></div><div  class = 'S7'><span>    Zero out the proton stoichiometr</span><span>ic coefficients from the temporary model for testing</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >modelH.S(protonMetBool,:)=0;</span></span></div></div></div><div  class = 'S7'><span>Check for duplicate columns, </span><span>upto protons, by detecting the columns of the  S matrix that are identical upto scalar multiplication.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >dupDetectMethod=</span><span style="color: rgb(170, 4, 249);">'FR'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >removeFlag=0;</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S9'><span style="white-space: normal"><span >[modelOut,removedRxnInd, keptRxnInd] = checkDuplicateRxn(modelH,dupDetectMethod,removeFlag,printLevel-1);</span></span></div><div  class = 'S10'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="F59D84D1" data-testid="output_13" data-width="420" data-height="73" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Checking for reaction duplicates by stoichiometry (up to orientation) ...
     Keep: 	BTNt2	btn[e] 	&lt;=&gt;	btn[c] 
Duplicate: 	BTNt4i	btn[e] 	-&gt;	btn[c] 
     Keep: 	GLCt1r	glc_D[e] 	&lt;=&gt;	glc_D[c] 
Duplicate: 	GLCt2_2	glc_D[e] 	&lt;=&gt;	glc_D[c] </div></div><div class="inlineElement eoOutputWrapper embeddedOutputsWarningElement" uid="0CCCAC2B" data-testid="output_14" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="diagnosticMessage-wrapper diagnosticMessage-warningType" style="white-space: normal; font-style: normal; color: rgb(255, 100, 0); font-size: 12px;"><div class="diagnosticMessage-messagePart" style="white-space: normal; font-style: normal; color: rgb(255, 100, 0); font-size: 12px;">Warning: Htg has more than one replicate</div><div class="diagnosticMessage-stackPart" style="white-space: normal; font-style: normal; color: rgb(255, 100, 0); font-size: 12px;"></div></div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="0B63A6FA" data-testid="output_15" data-width="420" data-height="143" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">     Keep: 	Htg		&lt;=&gt;	
Duplicate: 	HMR_1095		&lt;=&gt;	
     Keep: 	NACUP	nac[e] 	-&gt;	nac[c] 
Duplicate: 	NACDe	nac[c] 	-&gt;	nac[e] 
     Keep: 	ORNt4m	orn[m] + citr_L[c] 	&lt;=&gt;	orn[c] + citr_L[m] 
Duplicate: 	r0947	orn[m] + citr_L[c] 	-&gt;	orn[c] + citr_L[m] 
     Keep: 	THYMDt1	thymd[e] 	-&gt;	thymd[c] 
Duplicate: 	THYMDtr2	thymd[c] 	&lt;=&gt;	thymd[e] 
     Keep: 	VITD3tm	vitd3[m] 	-&gt;	vitd3[c] 
Duplicate: 	VITD3tm3	vitd3[c] 	-&gt;	vitd3[m] </div></div></div></div></div><div  class = 'S7'><span>Remove any du</span><span>plicate reactions from the stoichiometric matrix, but do not remove the protons.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">if </span><span >length(removedRxnInd)&gt;0</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    irrevFlag=0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    metFlag=0;</span><span style="color: rgb(2, 128, 9);">%dont remove the protons</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    model = removeRxns(model,model.rxns(removedRxnInd),irrevFlag,metFlag);</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div></div><div  class = 'S7'><span>Display statistics of the removed reactions</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">if </span><span >printLevel&gt;0</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    [nMet0,nRxn0]=size(modelOrig.S);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    [nMet,nRxn]=size(model.S);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%6s\t%6s\n'</span><span >,</span><span style="color: rgb(170, 4, 249);">'#mets'</span><span >,</span><span style="color: rgb(170, 4, 249);">'#rxns'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%6u\t%6u\t%s\n'</span><span >,nMet0,nRxn0,</span><span style="color: rgb(170, 4, 249);">' totals.'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%6u\t%6u\t%s\n'</span><span >,nMet0-nMet,nRxn0-nRxn,</span><span style="color: rgb(170, 4, 249);">' duplicate reactions upto protons removed.'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%6u\t%6u\t%s\n'</span><span >,nMet,nRxn,</span><span style="color: rgb(170, 4, 249);">' remaining.'</span><span >)</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div><div  class = 'S10'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="22B60F73" data-testid="output_16" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"> #mets	 #rxns</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="ED814FFD" data-testid="output_17" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">  8399	 13290	 totals.</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="75E053E3" data-testid="output_18" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">     0	     7	 duplicate reactions upto protons removed.</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="1BA91294" data-testid="output_19" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">  8399	 13283	 remaining.</div></div></div></div><div class="inlineWrapper"><div  class = 'S12'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%model size</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >[nMet,nRxn]=size(model.S);</span></span></div></div></div><h2  class = 'S3'><span>Heuristically identify exc</span><span>hange reactions and metabolites exclusively involved in exchange reactions</span></h2><div  class = 'S2'><span>An </span><span>external reacti</span><span>on is one tha</span><span>t is heuristically ide</span><span>ntified by a </span><span>s</span><span>ingle stoichiometric coefficient in the corresponding column of </span><span> </span><span>S, or an (abbreviated) reaction name matching a pattern (e.g. prefix EX_) or an external subsystem assignment. Any remaining reaction is assumed to be an internal reaction. If a reaction is not external then it is denoted an internal reaction. External reactants are exclusively involved in exchange reactions, and internal reactants otherwise. The findSExRxnInd function finds the external reactions in the model which export or import mass from or to the model, e.g. Exchange reactions, Demand reactions, Sink reactions.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">if </span><span >~isfield(model,</span><span style="color: rgb(170, 4, 249);">'SIntMetBool'</span><span >)  ||  ~isfield(model,</span><span style="color: rgb(170, 4, 249);">'SIntRxnBool'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >     model = findSExRxnInd(model,[],printLevel-1);</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div></div><h2  class = 'S1'><span>EXPECTED RESULTS</span></h2><div  class = 'S2'><span>In </span><span>the returned model, model.SIntRxnBool, is a boolean of reactions heuristically though to be mass balanced, while model.SIntMetBool is a boolean of metabolites heuristically though to be involved in mass balanced reactions.</span></div><h2  class = 'S1'><span>CAUTION</span></h2><div  class = 'S2'><span>The aforementioned assignments of external and internal reactions and reactants is the result of a heuristic and might result in one or more errors, either due to misspecification or because the names of external reactions and external subsystems often vary between laboratories. </span></div><h2  class = 'S3'><span>Find the reactions that are flux inconsistent </span></h2><div  class = 'S2'><span>Ultimately we seek to identify the set of stoichiometrically consistent reactions that are also flux consistent, with no bounds on reaction rates. However, finiding the stoichiometrically consistent subset can be demanding for large models so first we identify the subset of reactions that are flux consistent and focus on them.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >modelOrig=model;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >model.lb(~model.SIntRxnBool)=-1000;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >model.ub(~model.SIntRxnBool)= 1000;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">if </span><span >1</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">if </span><span >~isfield(model,</span><span style="color: rgb(170, 4, 249);">'fluxConsistentMetBool'</span><span >) || ~isfield(model,</span><span style="color: rgb(170, 4, 249);">'fluxConsistentRxnBool'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        param.epsilon=1e-4;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        param.modeFlag=0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        param.method=</span><span style="color: rgb(170, 4, 249);">'null_fastcc'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(2, 128, 9);">%param.method='fastcc';</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        [fluxConsistentMetBool,fluxConsistentRxnBool,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            fluxInConsistentMetBool,fluxInConsistentRxnBool,model]</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            = findFluxConsistentSubset(model,param,printLevel);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(2, 128, 9);">% Remove reactions that are flux inconsistent</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">if </span><span >any(fluxInConsistentRxnBool)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        irrevFlag=0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        metFlag=1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        model = removeRxns(model,model.rxns(fluxInConsistentRxnBool),irrevFlag,metFlag);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        [nMet0,nRxn0]=size(modelOrig.S);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        [nMet,nRxn]=size(model.S);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">if </span><span >printLevel&gt;0</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            fprintf(</span><span style="color: rgb(170, 4, 249);">'%s\n'</span><span >,</span><span style="color: rgb(170, 4, 249);">'-------'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            fprintf(</span><span style="color: rgb(170, 4, 249);">'%6s\t%6s\n'</span><span >,</span><span style="color: rgb(170, 4, 249);">'#mets'</span><span >,</span><span style="color: rgb(170, 4, 249);">'#rxns'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            fprintf(</span><span style="color: rgb(170, 4, 249);">'%6u\t%6u\t%s\n'</span><span >,nMet0,nRxn0,</span><span style="color: rgb(170, 4, 249);">' totals.'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            fprintf(</span><span style="color: rgb(170, 4, 249);">'%6u\t%6u\t%s\n'</span><span >,nMet0-nMet,nRxn0-nRxn,</span><span style="color: rgb(170, 4, 249);">' flux inconsistent reactions removed.'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            fprintf(</span><span style="color: rgb(170, 4, 249);">'%6u\t%6u\t%s\n'</span><span >,nMet,nRxn,</span><span style="color: rgb(170, 4, 249);">' remaining.'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            fprintf(</span><span style="color: rgb(170, 4, 249);">'%s\n'</span><span >,</span><span style="color: rgb(170, 4, 249);">'-------'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span><span style="color: rgb(14, 0, 255);">if </span><span >printLevel&gt;1</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >                </span><span style="color: rgb(14, 0, 255);">for </span><span >n=1:nRxn0</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >                    </span><span style="color: rgb(14, 0, 255);">if </span><span >fluxInConsistentRxnBool(n)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >                        fprintf(</span><span style="color: rgb(170, 4, 249);">'%15s\t%-100s\n'</span><span >,modelOrig.rxns{n},modelOrig.rxnNames{n})</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >                    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >                </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(2, 128, 9);">%revise model size</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        [nMet,nRxn]=size(model.S);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(2, 128, 9);">%Recompute</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(2, 128, 9);">%Heuristically identify exchange reactions and metabolites exclusively involved in exchange reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(2, 128, 9);">%finds the reactions in the model which export/import from the model</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(2, 128, 9);">%boundary i.e. mass unbalanced reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(2, 128, 9);">%e.g. Exchange reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(2, 128, 9);">%     Demand reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(2, 128, 9);">%     Sink reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        model = findSExRxnInd(model,[],0);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">if </span><span >printLevel&gt;0</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            fprintf(</span><span style="color: rgb(170, 4, 249);">'%s\n'</span><span >,</span><span style="color: rgb(170, 4, 249);">'------end------'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div><div  class = 'S10'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="0AF29E18" data-testid="output_20" data-width="420" data-height="227" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"> 12156	Total reactions
  5966	Reversible reactions.
  6190	Irreversible reactions.
 11488	Flux consistent reactions, without flipping.
   292	Flux inconsistent irreversible reactions, without flipping.
   376	Flux inconsistent reactions, without flipping.
 11783	Flux consistent reactions.
    81	Flux inconsistent reversible reactions left to flip.
 11785	Flux consistent reactions.
    79	Flux inconsistent reversible reactions left to flip.
 11787	Flux consistent reactions.
    77	Flux inconsistent reversible reactions left to flip.
 11789	Flux consistent reactions.
    45	Flux inconsistent reversible reactions left to flip.
 11791	Flux consistent reactions.
    35	Flux inconsistent reversible reactions left to flip.</div></div></div></div></div><h2  class = 'S3'><span>Fin</span><span>d mas</span><span>s leaks or siphons within the heuristically internal part, without using the bounds given by the model</span></h2><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">if </span><span >1</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    modelBoundsFlag=0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    leakParams.epsilon=1e-4;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    leakParams.method=</span><span style="color: rgb(170, 4, 249);">'dc'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    leakParams.theta=0.5;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    [leakMetBool,leakRxnBool,siphonMetBool,siphonRxnBool,leakY,siphonY,statp,statn] =</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        findMassLeaksAndSiphons(model,model.SIntMetBool,model.SIntRxnBool,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        modelBoundsFlag,leakParams,printLevel);</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div></div><h2  class = 'S3'><span>Fin</span><span>d the maximal set of reactions that are stoichiometrically consistent</span></h2><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">if </span><span >~isfield(model,</span><span style="color: rgb(170, 4, 249);">'SConsistentMetBool'</span><span >) || ~isfield(model,</span><span style="color: rgb(170, 4, 249);">'SConsistentRxnBool'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">if </span><span >strcmp(model.modelID,</span><span style="color: rgb(170, 4, 249);">'HMRdatabase2_00'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        massBalanceCheck=0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">else</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        massBalanceCheck=1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">if </span><span >1</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        [SConsistentMetBool,SConsistentRxnBool,SInConsistentMetBool,SInConsistentRxnBool,unknownSConsistencyMetBool,unknownSConsistencyRxnBool,model]</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            =findStoichConsistentSubset(model,massBalanceCheck,printLevel);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">else</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(2, 128, 9);">%print out problematic reactions to file</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        resultsFileName=[resultsPath filesep model.modelID];</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        [SConsistentMetBool,SConsistentRxnBool,SInConsistentMetBool,SInConsistentRxnBool,unknownSConsistencyMetBool,unknownSConsistencyRxnBool,model]</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            =findStoichConsistentSubset(model,massBalanceCheck,printLevel,resultsFileName);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >rxnBool=model.SInConsistentRxnBool &amp; model.SIntRxnBool;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">if </span><span >any(rxnBool)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">if </span><span >printLevel&gt;0</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        fprintf(</span><span style="color: rgb(170, 4, 249);">'%s\n'</span><span >,</span><span style="color: rgb(170, 4, 249);">'Stoichiometrically inconsistent heuristically non-exchange reactions:'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">for </span><span >n=1:nRxn</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">if </span><span >rxnBool(n)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            fprintf(</span><span style="color: rgb(170, 4, 249);">'%20s\t%50s\t%s\n'</span><span >,model.rxns{n},model.rxnNames{n},model.subSystems{n})</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">if </span><span >printLevel&gt;0</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        fprintf(</span><span style="color: rgb(170, 4, 249);">'%s\n'</span><span >,</span><span style="color: rgb(170, 4, 249);">'--------------'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >rxnBool=model.unknownSConsistencyRxnBool &amp; model.SIntRxnBool;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">if </span><span >any(rxnBool)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">if </span><span >printLevel&gt;0</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        fprintf(</span><span style="color: rgb(170, 4, 249);">'%s\n'</span><span >,</span><span style="color: rgb(170, 4, 249);">'Unknown consistency heuristically non-exchange reactions:'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">for </span><span >n=1:nRxn</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">if </span><span >rxnBool(n)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            fprintf(</span><span style="color: rgb(170, 4, 249);">'%20s\t%50s\t%s\n'</span><span >,model.rxns{n},model.rxnNames{n},model.subSystems{n})</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">if </span><span >printLevel&gt;0</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        fprintf(</span><span style="color: rgb(170, 4, 249);">'%s\n'</span><span >,</span><span style="color: rgb(170, 4, 249);">'--------------'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div></div><h2  class = 'S3'><span>San</span><span>ity check of stoichiometric and flux consistency of model</span><span> with open external reactions</span></h2><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">if</span><span >  all(model.SIntMetBool &amp; model.SConsistentMetBool)</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            &amp;&amp; nnz(model.SIntRxnBool &amp; model.SConsistentRxnBool)==nnz(model.SIntRxnBool)</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            &amp;&amp; all(model.fluxConsistentMetBool)</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            &amp;&amp; all(model.fluxConsistentRxnBool)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        [nMet,nRxn]=size(model.S);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">if </span><span >printLevel&gt;1</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            fprintf(</span><span style="color: rgb(170, 4, 249);">'%6s\t%6s\n'</span><span >,</span><span style="color: rgb(170, 4, 249);">'#mets'</span><span >,</span><span style="color: rgb(170, 4, 249);">'#rxns'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            fprintf(</span><span style="color: rgb(170, 4, 249);">'%6u\t%6u\t%s\n'</span><span >,nMet,nRxn,</span><span style="color: rgb(170, 4, 249);">' totals.'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            fprintf(</span><span style="color: rgb(170, 4, 249);">'%6u\t%6u\t%s\n'</span><span >,nnz(~model.SIntMetBool),nnz(~model.SIntRxnBool),</span><span style="color: rgb(170, 4, 249);">' heuristically exchange.'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        checksPassed=0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(2, 128, 9);">%Check that all heuristically non-exchange reactions are also stoichiometrically consistent</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(2, 128, 9);">%exchange reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        model.EXRxnBool=strncmp(</span><span style="color: rgb(170, 4, 249);">'EX_'</span><span >, model.rxns, 3)==1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(2, 128, 9);">%demand reactions going out of model</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        model.DMRxnBool=strncmp(</span><span style="color: rgb(170, 4, 249);">'DM_'</span><span >, model.rxns, 3)==1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(2, 128, 9);">%sink reactions going into or out of model</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        model.SinkRxnBool=strncmp(</span><span style="color: rgb(170, 4, 249);">'sink_'</span><span >, model.rxns, 5)==1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(2, 128, 9);">%all heuristic non-exchanges, i.e., supposedly all external reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        bool=~(model.EXRxnBool | model.DMRxnBool | model.SinkRxnBool);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">if </span><span >nnz(bool &amp; model.SIntRxnBool &amp; model.SConsistentRxnBool)==nnz(model.SConsistentRxnBool)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            checksPassed=checksPassed+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span><span style="color: rgb(14, 0, 255);">if </span><span >printLevel&gt;1</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >                fprintf(</span><span style="color: rgb(170, 4, 249);">'%6u\t%6u\t%s\n'</span><span >,nnz(model.SIntMetBool),nnz(model.SIntRxnBool),</span><span style="color: rgb(170, 4, 249);">' All internally stoichiometrically consistent. (Check 1: minimum cardinality of conservation relaxation vector.)'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(2, 128, 9);">%Check for mass leaks or siphons in the stoichiometrically consistent part</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(2, 128, 9);">%There should be no leaks or siphons in the stiochiometrically consistent part</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        modelBoundsFlag=0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        leakParams.epsilon=1e-4;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        leakParams.eta = getCobraSolverParams(</span><span style="color: rgb(170, 4, 249);">'LP'</span><span >, </span><span style="color: rgb(170, 4, 249);">'feasTol'</span><span >)*100;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        leakParams.method=</span><span style="color: rgb(170, 4, 249);">'dc'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        [leakMetBool,leakRxnBool,siphonMetBool,siphonRxnBool,leakY,siphonY,statp,statn]</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            =findMassLeaksAndSiphons(model,model.SConsistentMetBool,model.SConsistentRxnBool,modelBoundsFlag,leakParams,printLevel);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">if </span><span >nnz(leakMetBool)==0 &amp;&amp; nnz(leakRxnBool)==0 &amp;&amp; nnz(siphonMetBool)==0 &amp;&amp; nnz(siphonRxnBool)==0</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            checksPassed=checksPassed+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span><span style="color: rgb(14, 0, 255);">if </span><span >printLevel&gt;1</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >                fprintf(</span><span style="color: rgb(170, 4, 249);">'%6u\t%6u\t%s\n'</span><span >,nnz(leakMetBool | siphonMetBool),nnz(leakRxnBool | siphonRxnBool),</span><span style="color: rgb(170, 4, 249);">' No internal leaks or siphons. (Check 2: leak/siphon tests.)'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(2, 128, 9);">%Check that the maximal conservation vector is nonzero for each the</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(2, 128, 9);">%internal stoichiometric matrix</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        maxCardinalityConsParams.epsilon=1e-4;</span><span style="color: rgb(2, 128, 9);">%1/epsilon is the largest mass considered, needed for numerical stability</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        maxCardinalityConsParams.method = </span><span style="color: rgb(170, 4, 249);">'quasiConcave'</span><span >;</span><span style="color: rgb(2, 128, 9);">%seems to work the best, but sometimes infeasible</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        maxCardinalityConsParams.theta = 0.5;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        maxCardinalityConsParams.eta=getCobraSolverParams(</span><span style="color: rgb(170, 4, 249);">'LP'</span><span >, </span><span style="color: rgb(170, 4, 249);">'feasTol'</span><span >)*100;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        [maxConservationMetBool,maxConservationRxnBool,solution]=maxCardinalityConservationVector(model.S(model.SConsistentMetBool,model.SConsistentRxnBool), maxCardinalityConsParams);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">if </span><span >nnz(maxConservationMetBool)==size(model.S,1) &amp;&amp; nnz(maxConservationRxnBool)==nnz(model.SIntRxnBool)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            checksPassed=checksPassed+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span><span style="color: rgb(14, 0, 255);">if </span><span >printLevel&gt;1</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >                fprintf(</span><span style="color: rgb(170, 4, 249);">'%6u\t%6u\t%s\n'</span><span >,nnz(maxConservationMetBool),nnz(maxConservationRxnBool),</span><span style="color: rgb(170, 4, 249);">' All internally stoichiometrically consistent. (Check 3: maximim cardinality conservation vector.)'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(2, 128, 9);">%Check that each of the reactions in the model (with open external reactions) is flux consistent</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        modelOpen=model;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        modelOpen.lb(~model.SIntRxnBool)=-1000;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        modelOpen.ub(~model.SIntRxnBool)= 1000;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        param.epsilon=1e-4;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        param.modeFlag=0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        param.method=</span><span style="color: rgb(170, 4, 249);">'null_fastcc'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        [fluxConsistentMetBool,fluxConsistentRxnBool,fluxInConsistentMetBool,fluxInConsistentRxnBool,modelOpen] = findFluxConsistentSubset(modelOpen,param,printLevel-2);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">if </span><span >nnz(fluxConsistentMetBool)==size(model.S,1) &amp;&amp; nnz(fluxConsistentRxnBool)==size(model.S,2)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            checksPassed=checksPassed+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span><span style="color: rgb(14, 0, 255);">if </span><span >printLevel&gt;1</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >                fprintf(</span><span style="color: rgb(170, 4, 249);">'%6u\t%6u\t%s\n'</span><span >,nnz(fluxConsistentMetBool),nnz(fluxConsistentRxnBool),</span><span style="color: rgb(170, 4, 249);">' All flux consistent. (Check 4: maximim cardinality constrained right nullspace.)'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">if </span><span >checksPassed==4</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span><span style="color: rgb(2, 128, 9);">%save the model with open exchanges as the default generic</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span><span style="color: rgb(2, 128, 9);">%model</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            model=modelOpen;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span><span style="color: rgb(14, 0, 255);">if </span><span >printLevel&gt;0</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >                fprintf(</span><span style="color: rgb(170, 4, 249);">'%s\n'</span><span >,</span><span style="color: rgb(170, 4, 249);">'Open external reactions is stoichiometrically and flux consistent. A flux balance model generated from a reconstruction. GREAT!!!!'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >        save([resultsFileName </span><span style="color: rgb(170, 4, 249);">'_consistent.mat'</span><span >],</span><span style="color: rgb(170, 4, 249);">'model'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div></div><h2  class = 'S1'><span>R</span><span>EFE</span><span>RENCES</span></h2><div  class = 'S2'><span>Gevorgyan, A., Poolman, M. G., Fell D., Detection of stoichiometric inconsistencies in biomolecular models. Bioinformatics, 24(19):2245–51, 2008.</span></div><div  class = 'S2'><span>Flem</span><span>ing, R.M.T., et al., Cardinality optimisation in constraint-based modelling: Application to Recon 3D (submitted), 2017.</span></div><div  class = 'S2'><span>Brunk, E. et al. Recon 3D: A resource enabling a three-dimensional view of gene variation in human metabolism. (submitted) 2017.</span></div>
<br>
<!-- 
##### SOURCE BEGIN #####
%% *Convert a reconstruction into a flux balance analysis model* 
%% *Author: Ronan Fleming, Ines Thiele, University of Luxembourg*
%% *Reviewers:* 
%% INTRODUCTION
% Even with quality control during the reconstruction process, it is not appropriate 
% to assume that any reconstruction can be converted directly into a model and 
% used to make predictions. A model must satisfy certain assumptions before it 
% can be used to make reliable predictions. Depending on the type of model model, 
% these assumptions will be different. Each assumption should be chemically or 
% biologically motivated and expressed in an unambiguous manner and preferably 
% both intuitively and mathematically. Flux balance analysis is a mathematical 
% method widely used for studying genome-scale biochemical network. Here one aims 
% to predict steady-state reaction fluxes, where there is a balance between production 
% and consumption of each molecular species that is not exchanged across the specified 
% boundary of a system. In this situation, one might obtain erroneous predictions 
% if the system boundary is incorrectly specified. If a reconstruction contains 
% one or more supposedly mass balanced reactions, but which are actually not mass 
% balanced, such reactions in a model can lead to inadvertent leakage of a metabolite 
% from the model, in violation of mass balance. Similarly, when generating a model 
% for flux balance analysis, it is important to ensure that the network is flux 
% consistent, that is, each reaction can carry a non-zero steady state flux. 
% 
% Given a reconstruction with $$\hat{m}$$ reactants involved in $$\hat{n}$$ 
% reactions, this tutorial demonstrates a method to identify and extract the largest 
% subset of the reconstruction whose internal reactions are both stoichoimetrically 
% and flux consistent and whose external reactions are flux consistent. This model 
% is then mathematically consistent with the basic requirements for generation 
% of predictions using flux balance analysis. The identification of the component 
% of the reconstruction that does not satisfy the aforementioned modelling conditions 
% is also useful for targeting reconstruction effort towards resolving stoichiometric 
% inconsistency or resolving flux inconsistency. The example used in this tutorial 
% illustrates the process of extracting a model consistent with flux balance analsis, 
% from a ReconX reconstruction.
%% PROCEDURE
%% Select reconstruction to convert into a model and enter parameters
% Load the ReconX reconstruction, and save the original reconstruction in the 
% workspace, unless it is already loaded into the workspace. 

clear model
if ~exist('modelOrig','var')
    %select your own model, or use Recon2.0model instead
    if 1
        filename='Recon3D_301.mat'
        load(filename);
        model=Recon3D;
    else
        filename='Recon2.0model.mat';
        if exist('Recon2.0model.mat','file')==2
            model = readCbModel(filename);
        end
    end
    model.csense(1:size(model.S,1),1)='E';
    modelOrig = model;
else
    model=modelOrig;
end
%% 
% Set the level of printing, zero for silent, higher for more output.

printLevel=2;
%% 
% Choose the directory to place the results

basePath='~/work/sbgCloud/';
%resultsPath=[basePath '/programReconstruction/projects/recon2models/results/reconXs/' model.modelID];
resultsPath=[basePath '/courses/2019_Leiden_COBRA/practicalsDemo/Day4/' model.modelID];
resultsFileName=[resultsPath filesep model.modelID];
%% 
% Create and enter the folder for the results if it does not already exist

if ~exist(resultsPath,'dir')
    mkdir(resultsPath)
end
cd(resultsPath)
%% 
% Optionally create a diary to save the output in case it is very long, this 
% makes it easier to search, especially when debugging the process during the 
% early stages.

if 0
    diary([resultsFileName '_diary.txt'])
end
%% Overview some of the key properties of the reconstruction
% Noting the initial size of the reconstruction is useful for comparisons later 
% with subsets derived according to mathematical specifications.

[nMet,nRxn]=size(model.S);
fprintf('%6s\t%6s\n','#mets','#rxns')
fprintf('%6u\t%6u\t%s\n',nMet,nRxn,' totals.')
%% 
% Make sure the stoichiometric matrix is stored in a sparse format as this accelerates 
% computations with large networks

model.S=sparse(model.S);
%% Check in case the reconstruction is a model that is already ready for flux balance analysis
% There is no need to run this live script any further if the reconstruction 
% already satisfies the conditions necessary for flux balance analysis. That is 
% if all internal reactants and reactions are stoichiometrically consistent, and 
% all reactions are flux consistent, then the reconstruction satisfies the criteria 
% to designate it a model ready for flux balance analysis.
% 
% SIntMetBool                     m x 1 Boolean of metabolites heuristically 
% though to be involved in mass balanced reactions.
% 
% SIntRxnBool                     n x 1 Boolean of reactions heuristically though 
% to be mass balanced.
% 
% SConsistentMetBool        m x 1 Boolean vector indicating consistent mets
% 
% SConsistentRxnBool        n x 1 Boolean vector indicating consistent rxns
% 
% fluxConsistentMetBool     m x 1 Boolean vector indicating flux consistent 
% mets
% 
% fluxConsistentRxnBool     n x 1 Boolean vector indicating flux consistent 
% rxns

if all(isfield(model,{'SIntMetBool','SIntRxnBool','SConsistentMetBool',...
        'SConsistentRxnBool','fluxConsistentMetBool','fluxConsistentRxnBool'}))
    if all(model.SIntMetBool & model.SConsistentMetBool)...
            && nnz(model.SIntRxnBool & model.SConsistentRxnBool)==nnz(model.SIntRxnBool)...
            && all(model.fluxConsistentMetBool)...
            && all(model.fluxConsistentRxnBool)
        fullyStoichAndFluxConsistent=1;
        fprintf('%s\n','Reconstruction is a model that is already ready for flux balance analysis')
    end
    return
else
    fullyStoichAndFluxConsistent=0;
    fprintf('%s\n','Reconstruction must be tested to check if it is ready for flux balance analysis')
end
%% Manually remove certain reactions from the reconstruction 
% Before attempting to algorithmically remove stoichiometrically or flux inconsistent 
% supposed internal reactions from a reconstruction to generate a model, there 
% is an option to review the content of the reconstruction and manually identify 
% reactions for removal. That is, there are two options:
% 
% A. Skip manual review of reconstruction content. Move to the next step.
% 
% B. Review the content of the reconstruction and omit any reactions that are 
% assumed to be stoichiometrically or flux inconsistent. With respect to stoichiometric 
% inconsistency, such reactions may be obviously mass imbalanced and not satisfy 
% the heuristic conditions for indentification as an exernal reaction. Alternatively, 
% such reactions may be identified by a previous pass through of this tutorial 
% as being of unknown stoichometric consistent (model.unknownSConsistencyRxnBool(j)==1), 
% after the largest stoichiometrically consistent subset of the network has been 
% is identified. This is an iterative process where multiple rounds of identification 
% of the largest stoichiometrically consistent set and manual curation of the 
% remainder that is of unknown stoichiometric consistency is necessary.

if strcmp(filename,'Recon3.0model')
    modelOrig=model;
    if 0
        if 1
            %Rename some of the biomass reactions to make them more obviously exchange
            %reactions
            model.rxns{strcmp(model.rxns,'biomass_reaction')}= 'EX_biomass_reaction';
            model.rxns{strcmp(model.rxns,'biomass_maintenance')}= 'EX_biomass_maintenance';
            model.rxns{strcmp(model.rxns,'biomass_maintenance_noTrTr')}= 'EX_biomass_maintenance_noTrTr';
            
            %ATP hydrolysis is not imbalanced like all the other demand reactions so
            %give it a different accronym ATPM = ATP Maintenance
            bool=strcmp('DM_atp_c_',model.rxns);
            model.rxns{bool}='ATPM';
        end
        [model,removeMetBool,removeRxnBool] = manuallyAdaptRecon3(model,printLevel);
    else
        [model,removeMetBool,removeRxnBool] = manuallyAdaptRecon3Ines(model,printLevel);
    end
    [nMet0,nRxn0]=size(modelOrig.S);
    [nMet,nRxn]=size(model.S);
    if nMet0==nMet && nRxn0==nRxn && printLevel>0
        fprintf('%s\n','REPLACE_WITH_DASH_DASH- Manually removing rows and columns of the stoichiometric matrixREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASH')
        fprintf('%6s\t%6s\n','#mets','#rxns')
        fprintf('%6u\t%6u\t%s\n',nMet0,nRxn0,' totals.')
        fprintf('%6u\t%6u\t%s\n',nMet0-nMet,nRxn0-nRxn,' manually removed.')
        fprintf('%6u\t%6u\t%s\n',nMet,nRxn,' remaining.')
    end
end
%% Remove any trivial rows and columns of the stoichiometric matrix
% Remove any zero rows or columns of the stoichiometric matrix

modelOrig=model;
model=removeTrivialStoichiometry(model);
[nMet0,nRxn0]=size(modelOrig.S);
[nMet,nRxn]=size(model.S);
if nMet0==nMet && nRxn0==nRxn && printLevel>0
    fprintf('%s\n','REPLACE_WITH_DASH_DASH-Checking for Remove any trivial rows and columns of the stoichiometric matrixREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASH')
    fprintf('%6s\t%6s\n','#mets','#rxns')
    fprintf('%6u\t%6u\t%s\n',nMet0,nRxn0,' totals.')
    fprintf('%6u\t%6u\t%s\n',nMet0-nMet,nRxn0-nRxn,' duplicates removed.')
    fprintf('%6u\t%6u\t%s\n',nMet,nRxn,' remaining.')
end
%% 
% Check for duplicate columns by detecting the columns of the  S matrix that 
% are identical upto scalar multiplication.

modelOrig=model;
dupDetectMethod='FR';
dupDetectMethod='S';
removeFlag=0;
[modelOut,removedRxnInd, keptRxnInd] = checkDuplicateRxn(model,dupDetectMethod,removeFlag,printLevel-2);
%% 
% Remove any duplicate reactions, and uniquely involved reactants, from the 
% stoichiometric matrix. 

if length(removedRxnInd)>0
    irrevFlag=0;
    metFlag=1;
    %set all reactions reversible that are duplicates
    model.lb(removedRxnInd)=-model.ub(removedRxnInd);
    %remove duplicates
    model = removeRxns(model,model.rxns(removedRxnInd),irrevFlag,metFlag);
end
%% 
% Display the statistics on the duplicate reactions, 

[nMet0,nRxn0]=size(modelOrig.S);
[nMet,nRxn]=size(model.S);
if nMet0==nMet && nRxn0==nRxn && printLevel>0
    fprintf('%s\n','REPLACE_WITH_DASH_DASH-Remove any duplicate reactionsREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASH')
    [nMet0,nRxn0]=size(modelOrig.S);
    [nMet,nRxn]=size(model.S);
    fprintf('%6s\t%6s\n','#mets','#rxns')
    fprintf('%6u\t%6u\t%s\n',nMet0,nRxn0,' totals.')
    fprintf('%6u\t%6u\t%s\n',nMet0-nMet,nRxn0-nRxn,' duplicates removed.')
    fprintf('%6u\t%6u\t%s\n',nMet,nRxn,' remaining.')
end
%% Remove any duplicate reactions upto protons
% Remove reactions reactions that differ only in the number of protons involved 
% as substrates or products. Also remove exclusively involved reactants.
% 
% Save a temporary model for testing, before making any changes.

modelH=model;
%% 
% Find the proton indicies in different compartments. A proton, with index i, 
% is asumed to be represented by an abbreviation within model.mets{i} like h[*], 
% where * denotes the compartment symbol.

nMetChars=zeros(length(modelH.mets),1);
for m=1:length(modelH.mets)
    nMetChars(m,1)=length(modelH.mets{m});
end
protonMetBool=strncmp(modelH.mets,'h',1) & nMetChars==length('h[*]');
if printLevel>2
    disp(modelH.mets(protonMetBool))
end
%% 
% Zero out the proton stoichiometric coefficients from the temporary model for 
% testing

modelH.S(protonMetBool,:)=0;
%% 
% Check for duplicate columns, upto protons, by detecting the columns of the  
% S matrix that are identical upto scalar multiplication.

dupDetectMethod='FR';
removeFlag=0;
[modelOut,removedRxnInd, keptRxnInd] = checkDuplicateRxn(modelH,dupDetectMethod,removeFlag,printLevel-1);
%% 
% Remove any duplicate reactions from the stoichiometric matrix, but do not 
% remove the protons.

if length(removedRxnInd)>0
    irrevFlag=0;
    metFlag=0;%dont remove the protons
    model = removeRxns(model,model.rxns(removedRxnInd),irrevFlag,metFlag);
end
%% 
% Display statistics of the removed reactions

if printLevel>0
    [nMet0,nRxn0]=size(modelOrig.S);
    [nMet,nRxn]=size(model.S);
    fprintf('%6s\t%6s\n','#mets','#rxns')
    fprintf('%6u\t%6u\t%s\n',nMet0,nRxn0,' totals.')
    fprintf('%6u\t%6u\t%s\n',nMet0-nMet,nRxn0-nRxn,' duplicate reactions upto protons removed.')
    fprintf('%6u\t%6u\t%s\n',nMet,nRxn,' remaining.')
end
%model size
[nMet,nRxn]=size(model.S);
%% Heuristically identify exchange reactions and metabolites exclusively involved in exchange reactions
% An external reaction is one that is heuristically identified by a single stoichiometric 
% coefficient in the corresponding column of  S, or an (abbreviated) reaction 
% name matching a pattern (e.g. prefix EX_) or an external subsystem assignment. 
% Any remaining reaction is assumed to be an internal reaction. If a reaction 
% is not external then it is denoted an internal reaction. External reactants 
% are exclusively involved in exchange reactions, and internal reactants otherwise. 
% The findSExRxnInd function finds the external reactions in the model which export 
% or import mass from or to the model, e.g. Exchange reactions, Demand reactions, 
% Sink reactions.

if ~isfield(model,'SIntMetBool')  ||  ~isfield(model,'SIntRxnBool')
     model = findSExRxnInd(model,[],printLevel-1);
end
%% EXPECTED RESULTS
% In the returned model, model.SIntRxnBool, is a boolean of reactions heuristically 
% though to be mass balanced, while model.SIntMetBool is a boolean of metabolites 
% heuristically though to be involved in mass balanced reactions.
%% CAUTION
% The aforementioned assignments of external and internal reactions and reactants 
% is the result of a heuristic and might result in one or more errors, either 
% due to misspecification or because the names of external reactions and external 
% subsystems often vary between laboratories. 
%% Find the reactions that are flux inconsistent 
% Ultimately we seek to identify the set of stoichiometrically consistent reactions 
% that are also flux consistent, with no bounds on reaction rates. However, finiding 
% the stoichiometrically consistent subset can be demanding for large models so 
% first we identify the subset of reactions that are flux consistent and focus 
% on them.

modelOrig=model;
model.lb(~model.SIntRxnBool)=-1000;
model.ub(~model.SIntRxnBool)= 1000;
if 1
    if ~isfield(model,'fluxConsistentMetBool') || ~isfield(model,'fluxConsistentRxnBool')
        param.epsilon=1e-4;
        param.modeFlag=0;
        param.method='null_fastcc';
        %param.method='fastcc';
        [fluxConsistentMetBool,fluxConsistentRxnBool,...
            fluxInConsistentMetBool,fluxInConsistentRxnBool,model]...
            = findFluxConsistentSubset(model,param,printLevel);
    end
    % Remove reactions that are flux inconsistent
    if any(fluxInConsistentRxnBool)
        irrevFlag=0;
        metFlag=1;
        model = removeRxns(model,model.rxns(fluxInConsistentRxnBool),irrevFlag,metFlag);
        [nMet0,nRxn0]=size(modelOrig.S);
        [nMet,nRxn]=size(model.S);
        
        if printLevel>0
            fprintf('%s\n','REPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASH-')
            fprintf('%6s\t%6s\n','#mets','#rxns')
            fprintf('%6u\t%6u\t%s\n',nMet0,nRxn0,' totals.')
            fprintf('%6u\t%6u\t%s\n',nMet0-nMet,nRxn0-nRxn,' flux inconsistent reactions removed.')
            fprintf('%6u\t%6u\t%s\n',nMet,nRxn,' remaining.')
            fprintf('%s\n','REPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASH-')
            if printLevel>1
                for n=1:nRxn0
                    if fluxInConsistentRxnBool(n)
                        fprintf('%15s\t%-100s\n',modelOrig.rxns{n},modelOrig.rxnNames{n})
                    end
                end
            end
        end
        %revise model size
        [nMet,nRxn]=size(model.S);
        
        %Recompute
        %Heuristically identify exchange reactions and metabolites exclusively involved in exchange reactions
        %finds the reactions in the model which export/import from the model
        %boundary i.e. mass unbalanced reactions
        %e.g. Exchange reactions
        %     Demand reactions
        %     Sink reactions
        
        model = findSExRxnInd(model,[],0);
        if printLevel>0
            fprintf('%s\n','REPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHendREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASH')
        end
    end
end
%% Find mass leaks or siphons within the heuristically internal part, without using the bounds given by the model

if 1
    modelBoundsFlag=0;
    leakParams.epsilon=1e-4;
    leakParams.method='dc';
    leakParams.theta=0.5;
    [leakMetBool,leakRxnBool,siphonMetBool,siphonRxnBool,leakY,siphonY,statp,statn] =...
        findMassLeaksAndSiphons(model,model.SIntMetBool,model.SIntRxnBool,...
        modelBoundsFlag,leakParams,printLevel);
end
%% Find the maximal set of reactions that are stoichiometrically consistent

if ~isfield(model,'SConsistentMetBool') || ~isfield(model,'SConsistentRxnBool')
    if strcmp(model.modelID,'HMRdatabase2_00')
        massBalanceCheck=0;
    else
        massBalanceCheck=1;
    end
    if 1
        [SConsistentMetBool,SConsistentRxnBool,SInConsistentMetBool,SInConsistentRxnBool,unknownSConsistencyMetBool,unknownSConsistencyRxnBool,model]...
            =findStoichConsistentSubset(model,massBalanceCheck,printLevel);
    else
        %print out problematic reactions to file
        resultsFileName=[resultsPath filesep model.modelID];
        [SConsistentMetBool,SConsistentRxnBool,SInConsistentMetBool,SInConsistentRxnBool,unknownSConsistencyMetBool,unknownSConsistencyRxnBool,model]...
            =findStoichConsistentSubset(model,massBalanceCheck,printLevel,resultsFileName);
    end
end
    
rxnBool=model.SInConsistentRxnBool & model.SIntRxnBool;
if any(rxnBool)
    if printLevel>0
        fprintf('%s\n','Stoichiometrically inconsistent heuristically non-exchange reactions:')
    end
    for n=1:nRxn
        if rxnBool(n)
            fprintf('%20s\t%50s\t%s\n',model.rxns{n},model.rxnNames{n},model.subSystems{n})
        end
    end
    if printLevel>0
        fprintf('%s\n','REPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASH')
    end
end
    
rxnBool=model.unknownSConsistencyRxnBool & model.SIntRxnBool;
if any(rxnBool)
    if printLevel>0
        fprintf('%s\n','Unknown consistency heuristically non-exchange reactions:')
    end
    for n=1:nRxn
        if rxnBool(n)
            fprintf('%20s\t%50s\t%s\n',model.rxns{n},model.rxnNames{n},model.subSystems{n})
        end
    end
    if printLevel>0
        fprintf('%s\n','REPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASH')
    end
end
%% Sanity check of stoichiometric and flux consistency of model with open external reactions

    if  all(model.SIntMetBool & model.SConsistentMetBool)...
            && nnz(model.SIntRxnBool & model.SConsistentRxnBool)==nnz(model.SIntRxnBool)...
            && all(model.fluxConsistentMetBool)...
            && all(model.fluxConsistentRxnBool)
        
        [nMet,nRxn]=size(model.S);
        if printLevel>1
            fprintf('%6s\t%6s\n','#mets','#rxns')
            fprintf('%6u\t%6u\t%s\n',nMet,nRxn,' totals.')
            fprintf('%6u\t%6u\t%s\n',nnz(~model.SIntMetBool),nnz(~model.SIntRxnBool),' heuristically exchange.')
        end
        
        checksPassed=0;
        %Check that all heuristically non-exchange reactions are also stoichiometrically consistent
        
        %exchange reactions
        model.EXRxnBool=strncmp('EX_', model.rxns, 3)==1;
        %demand reactions going out of model
        model.DMRxnBool=strncmp('DM_', model.rxns, 3)==1;
        %sink reactions going into or out of model
        model.SinkRxnBool=strncmp('sink_', model.rxns, 5)==1;
        %all heuristic non-exchanges, i.e., supposedly all external reactions
        bool=~(model.EXRxnBool | model.DMRxnBool | model.SinkRxnBool);
        if nnz(bool & model.SIntRxnBool & model.SConsistentRxnBool)==nnz(model.SConsistentRxnBool)
            checksPassed=checksPassed+1;
            if printLevel>1
                fprintf('%6u\t%6u\t%s\n',nnz(model.SIntMetBool),nnz(model.SIntRxnBool),' All internally stoichiometrically consistent. (Check 1: minimum cardinality of conservation relaxation vector.)');
            end
        end
        
        %Check for mass leaks or siphons in the stoichiometrically consistent part
        %There should be no leaks or siphons in the stiochiometrically consistent part
        modelBoundsFlag=0;
        leakParams.epsilon=1e-4;
        leakParams.eta = getCobraSolverParams('LP', 'feasTol')*100;
        leakParams.method='dc';
        [leakMetBool,leakRxnBool,siphonMetBool,siphonRxnBool,leakY,siphonY,statp,statn]...
            =findMassLeaksAndSiphons(model,model.SConsistentMetBool,model.SConsistentRxnBool,modelBoundsFlag,leakParams,printLevel);
        
        if nnz(leakMetBool)==0 && nnz(leakRxnBool)==0 && nnz(siphonMetBool)==0 && nnz(siphonRxnBool)==0
            checksPassed=checksPassed+1;
            if printLevel>1
                fprintf('%6u\t%6u\t%s\n',nnz(leakMetBool | siphonMetBool),nnz(leakRxnBool | siphonRxnBool),' No internal leaks or siphons. (Check 2: leak/siphon tests.)');
            end
        end
        
        %Check that the maximal conservation vector is nonzero for each the
        %internal stoichiometric matrix
        maxCardinalityConsParams.epsilon=1e-4;%1/epsilon is the largest mass considered, needed for numerical stability
        maxCardinalityConsParams.method = 'quasiConcave';%seems to work the best, but sometimes infeasible
        maxCardinalityConsParams.theta = 0.5;
        maxCardinalityConsParams.eta=getCobraSolverParams('LP', 'feasTol')*100;
        [maxConservationMetBool,maxConservationRxnBool,solution]=maxCardinalityConservationVector(model.S(model.SConsistentMetBool,model.SConsistentRxnBool), maxCardinalityConsParams);
        
        if nnz(maxConservationMetBool)==size(model.S,1) && nnz(maxConservationRxnBool)==nnz(model.SIntRxnBool)
            checksPassed=checksPassed+1;
            if printLevel>1
                fprintf('%6u\t%6u\t%s\n',nnz(maxConservationMetBool),nnz(maxConservationRxnBool),' All internally stoichiometrically consistent. (Check 3: maximim cardinality conservation vector.)');
            end
        end
        
        %Check that each of the reactions in the model (with open external reactions) is flux consistent
        modelOpen=model;
        modelOpen.lb(~model.SIntRxnBool)=-1000;
        modelOpen.ub(~model.SIntRxnBool)= 1000;
        param.epsilon=1e-4;
        param.modeFlag=0;
        param.method='null_fastcc';
        [fluxConsistentMetBool,fluxConsistentRxnBool,fluxInConsistentMetBool,fluxInConsistentRxnBool,modelOpen] = findFluxConsistentSubset(modelOpen,param,printLevel-2);
        
        if nnz(fluxConsistentMetBool)==size(model.S,1) && nnz(fluxConsistentRxnBool)==size(model.S,2)
            checksPassed=checksPassed+1;
            if printLevel>1
                fprintf('%6u\t%6u\t%s\n',nnz(fluxConsistentMetBool),nnz(fluxConsistentRxnBool),' All flux consistent. (Check 4: maximim cardinality constrained right nullspace.)');
            end
        end
        
        if checksPassed==4
            %save the model with open exchanges as the default generic
            %model
            model=modelOpen;
            if printLevel>0
                fprintf('%s\n','Open external reactions is stoichiometrically and flux consistent. A flux balance model generated from a reconstruction. GREAT!!!!');
            end
        end
        save([resultsFileName '_consistent.mat'],'model')
    end
%% REFERENCES
% Gevorgyan, A., Poolman, M. G., Fell D., Detection of stoichiometric inconsistencies 
% in biomolecular models. Bioinformatics, 24(19):2245–51, 2008.
% 
% Fleming, R.M.T., et al., Cardinality optimisation in constraint-based modelling: 
% Application to Recon 3D (submitted), 2017.
% 
% Brunk, E. et al. Recon 3D: A resource enabling a three-dimensional view of 
% gene variation in human metabolism. (submitted) 2017.
##### SOURCE END #####
-->
</div></body></html>